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We present a systematic investigation of attractive binary mixtures in presence of both spin- and mass- 
imbalance in one dimensional setups described by the Hubbard model. After discussing typical cold atomic 
experimental realizations and the relation between microscopic and effective parameters, we study several 
many-body features of trapped Fermi-Fermi and Bose-Bose mixtures such as density profiles, momentum dis- 
tributions and correlation functions by means of numerical density-matrix-renormalization-group and Quantum 
Monte Carlo simulations. In particular, we focus on the stability of Fulde-Ferrell-Larkin-Ovchinnikov, dimer 
and trimer fluids in inhomogeneous situations, as typically realized in cold gas experiments due to the harmonic 
confinement. We finally consider possible experimental signatures of these phases both in the presence of a 
finite polarization and of a finite temperature. 

PACS numbers: 71.10.Pm, 05.30.Jp, 37.10.Jk 



I. INTRODUCTION 

The study of superfluidity of either bosons or fermions has 
been a central topic in the field of ultracold atomic gases, start- 
ing from seminal experimental studies on the Bose-Einstein 
condensation of bosons m^ and continuing with investiga- 
tions of the BCS-BEC crossover |0-[2l]- More recent exper- 
iments with two-component Fermi gases have addressed the 
case of a finite population imbalanced, both in three spatial 
dimensions Il8l 4llll and in one dimension III2I1 . Among the 
goals of these experiments is the search for the transition from 
a fully paired superfluid phase to the normal state and, in par- 
ticular, for competing pairing states that survive a finite po- 
larization. These include the Sarma state J|13 |or the Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO) state U3> ll^ll • 

With the advances optical lattice engineering III6I1 . it has 
also become possible to investigate fermionic pairing states 
within the framework of the so-called Hubbard model. For 
repulsive onsite interactions the hope is to reach sufficiently 
low temperatures to search for pairing states away from half 
filling III7I419I] , whereas in the case of attractive interactions, 
there is a natural tendency to pair formation. Besides ex- 
periments with homonuclear mixtures, where the pseudo-spin 
degree of freedom arises from preparing atoms in different 
hyper-fine states, there is also the possibility of working with 
heteronuclear mixtures such as, for instance, the ^^K-^Li sys- 
tem |20-[23n. In that case, one deals both finite mass imbal- 
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ance and population imbalance. On the other hand, bosonic 
gases in optical lattices have provided the first example of a 
quantum phase transition from a Mott-insulator to a superfluid 
state in the context of ultra cold atoms ll24ll . which was subse- 
quently observed in one dimensional setups ||25l |2fl] as well. 
In analogy with the fermionic case, both homonuclear and 
heteronuclear bosonic mixtures have been realized [I27l - l29ll ; 
from a theoretical viewpoint, the two-species Bose-Hubbard 
model encompasses a remarkably rich physics, ranging from 
super-counterflow and antiferromagnetic phases in case of in- 
terspecies repulsion [30-32] to pair superfluidity and density- 
wave instabilities in the attractive regime [|33| - |35|] . 

From a theoretical point of view, it is of advantage to con- 
sider the one-dimensional case, for which both powerful an- 
alytical and numerical methods are available that can pro- 
vide us with practically exact answers in many regimes. In 
ultra-cold atomic gases, it is also perfectly possible to real- 
ize one-dimensional systems experimentally, both for bosons 
1I25I] and fermions (ij], adding strong motivation to study this 
case. Moreover one-dimensional systems have been demon- 
strated to harbor unconventional superfluid phases which are 
the central topic of this work. In our work, we are interested in 
three types of superfluid states of two-component systems and 
their existence in the attractive, asymmetric Hubbard model: 
(i) the more conventional fully paired phase, which can be 
considered as a fluid of Cooper pairs, (ii) superfluids of larger 
composite objects such as trimers, and (iii) the FFLO state 
in which a superfluid of pairs with finite center-of-mass mo- 
mentum appears. The first two cases, i.e., the equal density 
(or dimer fluid) and the trimer fluid, can be realized in both 
Fermi-Fermi mixtures ||36l |37ll and Bose-Bose mixtures in 



one-dimension 1331] while the FFLO state exists only in the 
case of two-component Fermi mixtures (see [.3811 and refer- 
ences therein). In our study we use numerically exact meth- 
ods - density matrix renormalization group (DMRG) [39^-41] 
and stochastic series expansion Quantum Monte Carlo [421 - 
to probe the stability of the above cited phases at a finite popu- 
lation and mass imbalance in the presence of a harmonic trap- 
ping potential. Our study complements previous works that 
mostly focussed on the homogeneous case ll34l - [37l l43l - [45ll . 
Moreover, we incorporate aspects that are typical of many ex- 
periments, namely the possibility of tuning the ratio of the ef- 
fective masses through the depth of the optical lattice; the fact 
that in a mass imbalanced system the two components usually 
experience a different trapping potential; the difficulty in con- 
trolling the population of the two components to an arbitrary 
degree of precision, resulting typically in a non-zero global 
polarization. Finally, we also study the effect of a small but 
finite temperature on the shell structure and coherence prop- 
erties of trapped mass-imbalanced two-component gases. 

Here we summarize our main results. In the case of 
fermions and for the bulk system it is well-known that the fully 
paired phase is the ground state of a two-component Fermi gas 
away from half filling ll46ll . In the presence of a harmonic trap 
and at a finite population imbalance, the fully paired phase 
can only survive in the wings of the particle cloud at small 
polarizations ||47| - |5Q|i . as recently observed in Ref. 12. Here 
we show that adding mass imbalance to the system stabilizes 
the fully paired phase in the trap: if the light particles are the 
majority species, the fully paired phase occupies the center 
of the trapped system, in agreement with the grand-canonical 
phase diagram of the asymmetric Hubbard model at a finite 
spin polarization ||3fl |4^ |45|] . We further show that a small 
temperature does not destroy the equal density phase. The 
main conclusion is therefore that it is not necessary to en- 
force the condition of a perfectly balanced gas to observe and 
study properties of a dimer fluid. This result applies to both 
fermions and hard-core bosons. 

A particular feature of mass-imbalanced systems in one di- 
mension is that they allow for the formation of bound states 
of more than two components, which in the two- species Hub- 
bard model with equal masses is forbidden. This can happen 
for both fermions ||36l llvl |45l] and bosons [37]. The existence 
of these states requires commensurate ratios of the densities 
of heavy particles (n^) versus light ones (n^). For instance, 
a trimer fluid can only be stable if 2n^ = n^ |13^[37|], where 
Ua- = Nfj/L is the density of the pseudo-spin a =t, i com- 
ponent in a system of length L with A^^r particles of species cr. 
For the homogeneous system, this sets the global polarization 
to 



(iVt-7V;)/Ar = -l/3 



(1) 



(A' = 7V^ + A^;) and it is in fact possible to stabilize this phase 
in a harmonic trap by enforcing this condition globally as was 
shown in Ref. Hat sufficiently low densities. Here we investi- 
gate the fate of the trimer fluid upon deviating from the global 
polarization of P = —1/3, both for the case of fermions and 
bosons, finding that the trimer fluid rapidly gets pushed into 
the wings of the gas and subsequently disappears even for very 



small deviations from Eq. [T] Therefore, one needs to fine- 
tune the global polarization to precisely P = —1/3 to ensure 
that a large fraction of particles participates in this state. A 
difference in the trapping potential for the two components, 
however, has little effect on the trimer phase. For the case of 
bosons, we demonstrate that the formation of the trimer fluid 
leaves clear fingerprints in the momentum distribution func- 
tion, a quantity that is easily accessible in experiments. 

Finally, concerning the FFLO state, we show that it sur- 
vives the addition of the harmonic trap and the mass imbal- 
ance. Similar to the case of a system with mass imbalance 
only ll47l - [5in . this phase occupies the center of the trap in 
a large parameter regime, corresponding in particular to the 
case in which heavy particles are the majority ones. These 
numerically exact observations are consistent with the known 
phase diagrams for ho mog eneous system via a local-density 
approximation iMHIUD. 

The plan of the paper is the following. First, in Sec. |II]we 
discuss typical experimental settings and conditions for two 
examples: (i) a heteronuclear Fermi-Fermi mixture such as 
^^K and ^Li and (ii) a homonuclear Bose-Bose mixture, where 
the mass-imbalance arises due to a spin-dependent optical lat- 
tice. This discussion serves to guide our numerical study. In 
Section [nil we provide an extensive analysis of the density 
profiles of a mass-imbalanced Fermi mixture and we present 
a state diagram for the various shell structures. In Sec. |IV] we 
study both hardcore and softcore bosons and probe the stabil- 
ity of the dimer and trimer fluid. We conclude each of these 
section with a detailed summary of the results, while the main 
aspects are summarized again in Sec.lVl 



II. EXPERIMENTAL SET-UPS FOR MASS-IMBALANCED 
SYSTEMS IN OPTICAL LATTICES 

The purpose of this section is to describe the experimen- 
tal set-ups to realize mass-imbalanced systems described by 
the asymmetric Hubbard model. In this context, we dis- 
cuss heteronuclear mixtures, focussing on the example of the 
Fermi-Fermi system ^^K and ^Li. The case of a homonu- 
clear Fermi-Fermi mixture for which the mass-imbalance is 
induced by spin-dependent optical lattices has been described 
in 152]. Then we turn to homonuclear Bose mixtures in a spin- 
dependent optical lattice. 



A. Model Hamiltonian 

Our numerical analysis will be based on the attractive, 
asymmetric ID Hubbard model: 



L-l 



l^l,a 1=1 

L 

+ ^ V^{l-L/2)\i, (2) 



l=l,a 



where cj^ creates a fermion with a pseudo-spin g =t, i at 

site /, UicT = cl^Cia^ ^i = ^i'\ + ^li is the local density, t^j 
is the hopping parameter with an explicit dependence on the 
pseudo-spin index, and C/ < is the attractive onsite inter- 
action energy. We set the lattice spacing to unity and impose 
open boundary conditions. We add a harmonic confining po- 
tential parameterized by constants Va- • 

Without loss of generality we will consider t^ > t^. A pos- 
itive (negative) polarization P corresponds to a majority of 
light (heavy) particles. We also introduce the local magneti- 
zation through 

{5f> = K^-n,^>/2. (3) 

We further define the ratio of the trapping potentials as 

n = n/n • (4) 

B. Heteronuclear Fermi-Fermi mixtures 
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FIG. 1: An array of independent one-dimensional optical dipole 
traps is created by intersecting two perpendicular standing waves 
(red beams). A third orthogonal optical standing wave (blue beams) 
generates a repulsive lattice potential along the traps. The density 
profiles of trapped atoms can be observed by an imaging setup from 
a transverse direction. 



The proposed experimental setup is schematically shown in 
Fig. [U Two pairs of counterpropagating laser beams with or- 
thogonal linear polarization are orthogonally intersected. This 
creates a 2D array of individual atom traps in the shape of one- 
dimensional tubes. In this example we choose the wavelength 
to be A*^^P = 1024 nm, far red detuned from the atomic tran- 
sitions. Considering a two- species fermionic mixture of ^Li 
and ^^K, we take into account the transition wavelengths for 
the D-lines (671.0 nm for lithium and 770.1 nm / 767.7 nm for 
potassium). In order to realize one-dimensional optical lat- 
tices we superimpose a pair counterpropagating laser beams 
perpendicular trapping beams and blue detuned with respect 
to the wavelengths of the D lines. In the following subsec- 
tions, we will discuss the resulting trapping potentials Va, the 
hopping parameters to- and the onsite interaction U. 



1. Trapping potentials 

The ratio of the optical trapping potentials y^^^^P {i =Li,K) 
is a fixed number that amounts for the given wavelengths to 



F^rV^x"'' - 0-44. 



(5) 



Note that V^ = 1/l[ and V^ = F^ • Due to the difference 
in the optical potentials and the different masses the respective 
trapping frequencies are generally different: 
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(6) 



where wi is the waist of the trapping beams and Vq\^^ is the 
trap depth in the center of the crossed beam setup. 

In general the wavelength for the trap can be chosen such 
that the difference in the optical potentials for the two atomic 
species compensates the mass difference such that equal trap- 
ping frequencies are obtained JSm. This so called magic 
wavelength, A^^^, occurs for the mixture of ^Li- ^^K at 
799.9 nm. 
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FIG. 2: The normalized lattice depth: At the magic wavelength 






661.3 nm) the lattice depths are equal for Li and K (solid 



line). Two nanometer above (dotted) and below (dashed) the lattice 
depths differ, while the ratio is inverted at the magic wavelength. 
Using a beam waist of 200 /xm less than one Watt of laser power is 
sufficient to achieve relevant lattice depths. 



2. The optical lattice 



The one-dimensional optical lattice is realized by superim- 
posing a pair of counterpropagating laser beams perpendicular 
to the beams that form the trapping array. For a given inten- 
sity of the lasers creating the optical lattice, potassium and 
lithium will experience optical lattice potentials of a different 
depths as illustrated in Fig.[2l The experimentally controllable 
parameters are the amplitudes of the standing wave forming 
the optical lattice, i.e., the lattice depths V^^f. These, and the 
masses determine the respective frequencies: 



(7) 



It is useful to introduce the atomic recoil energies as 



En^i 



2mi 27X110? 



(8) 



where k = 27r/A^^ is the wave vector of the standing wave, 
A^^ is the laser wavelength, and a = A^^/2 is the lattice spac- 
ing. This allows us to express the frequencies as: 



^i = 






(9) 



For wavelengths of the lattice laser blue detuned to the atomic 
transitions there exists a single magic wavelength, A^, at 
which the lattice depths Vq^^ normalized to the respective re- 
coil energies ER^i are equal: 



ysw 



= Fo, 



For the ^Li-^^K mixture this occurs at 661.3 nm. In Table IJ 
we summarize the magic wavelengths and the ratios of the 
optical trapping potentials at A*^^p = 1024 nm for different 
combinations of alkaline atomic species. As is obvious from 
the table, in all cases, the trapping potentials differ. 

We assume a separable three-dimensional lattice potential 
of the form 

V{x,y,z) = V\\{x)^V^{y)^V^{z) (10) 

= Fo,|| sin^(/cx) + Vo^± [sm'^{ky) + sin^ {kz)] . 

This allows us to calculate the hopping matrix elements from 
the one-dimensional Mathieu equation. The result is shown 
in Fig. [3] (solid lines). In the limit of a deep lattice, one can 
obtain an analytical expression llM : 
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V ER,^ 

(11) 

with i = K^Li and recoil energies ER^i. This result is also 
included in Fig. [3] (dashed lines), and it agrees with the exact 
solution for large Vq"^ or respectively laser power Pgw 

The ratio of the tunneling energies of ^Li and ^^K is shown 
in Fig.m While the bare mass ratio is mLi/rriK ~ 0.15, we 
deduce that the ratio of the actual hopping matrix elements 
further depends on the laser power Psw through the lattice 
depth Vq7. This is illustrated in Fig. [4] Thus in general, the 
ratio t^/t^ can be tuned and is not solely given by the bare 
mass ratio. In our numerical analysis of a two-component 
Fermi mixture, we will work at t^/t^ = 0.3, corresponding to 
thi/tK = 3. This requires to detune the wavelength slightly 
above the magic wavelength A^ (compare Fig.JH) 

The onsite repulsion U between two atoms of different 
species is given by Il54ll : 



U = ^^^^^ I dr'Mx,y,z)f\^,{x,y,z)f (12) 
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FIG. 3: The tunneling energies at the magic wavelength are the same 
for ^Li and ^°K (solid line). The asymptotic expression (dashed line) 
only describes the exact solution of the Mathieu equation at suffi- 
ciently high lattice depths (see Eq. ilTl). The maximum at low lattice 
depths is not physical. 
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FIG. 4: Ratio of the tunneling energies for ^Li and ^°K for different 
lattice wavelengths. At the magic wavelength the ratio is constant at 
the inverse mass ratio (solid line). Two nanometers above and below 
the magic wave-length the ratio of the tunneling energies differs from 
the bare mass ratio (dashed lines). By increasing the lattice power 
the ratio of the tunneling energies can be tuned over a wide range 
that can be set by the choice of the wavelength. The DMRG and 
SSES simulations presented in this work are done for a ratio of the 
tunneling energies thi/tK = 3|. 



In the expression, ijja- are the Wannier functions , jiij is the 
reduced mass, and a^^j is the scattering length that could be 
different for two atomic species, hyperfine flavors, or Bose- 
Fermi mixtures. 

For the evaluation of U we assume a deep lattice Vq^ ^ 
Er^i and that the particles occupy the lowest band. We can 
then approximate the Wannier functions by the ground state 
functions of a ID harmonic oscillator: 



V^i,||(^) 



(^a,2||)-V^exp(-xV2a,^l|) 

(^a|^)-i/^exp(-^V2a^^) 

(^a^^)-V4exp(-zV2a^^) 



(13) 



TABLE I: Magic wavelengths resulting in identical trapping frequen- 
cies or lattice potentials for different combinations of alkaline atoms. 
The calculation takes both D lines for each species into account. 





T/trap /T/trap 

at 1064 nm 


(nm) 


\sw 

(nm) 


6Li40K 


0.44 


799.9 


661.3 


^Li^^Rb 


0.4 


804.8 


666.3 


^Li^^a^s 


0.23 


905.1 


666.7 


23Na"°K 


0.39 


- 


481.1 


^^Na^^Rb 


0.35 


945.8 


551.7 


40^87^,^ 


0.9 


807.3 


655.4 
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Notice that for soft-core bosons, one needs to consider excited 
states. 

A simple and straightforward Gaussian integration yields: 
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This expression can be simplified if we work at the magic 
wavelength for the lattice, since for this case we have: 



1/4- 
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Expressing the energies in units of one of the recoil energies, 
we obtain 
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In the extreme ID limit, this entire expression can be sim- 
plified even further, since the Wannier functions become Dirac 
deltas in the transverse direction: 



C. Homonuclear two-component Bose gases in spin-dependent 
optical lattices 



In the case of bosons the Hubbard Hamiltonian contains as 
well intraspecies interactions: 

L 



i=l 



^t 



Here 6]^ creates a boson of type g =t,i on site i. In this 
case, mass imbalance can be realized either in heteronuclear 
mixtures (such as ^^Rb/^^K 11271 ISSH ) or in hyperfine mixtures 
in spin-dependent optical lattices (such as in the case of ^^Rb 
as realized in a number of recent experiments [|28l |29| . We 
will focus the following short discussion on the latter example, 
in which the Hubbard parameters can be tuned by the lattice 
depth, by the optical lattice wavelength, and by an interspecies 
Eeshbach resonance. We consider ^^Rb in the hyperfine mix- 
ture of I t) = |F = 1, rriF = 1) and | |) = |F = 1, m^ = 0) 
hyperfine states, which are collisionally stable and for which 
several Eeshbach resonances have been identified 11561] . State- 
dependent optical lattices can be easily realized when the op- 
tical lattice wavelength is close to the Dl and D2 lines - 
795 and 780 nm respectively IdtIi . We consider a highly 
anisotropic optical lattice, with wavelength A^ = 830 nm and 
depth Vo^± = 40 £^^ for the transverse components defining 
ID tubes, and wavelength A|| = 784 nm and variable depth 
Fo,|| for the longitudinal component along the tubes. Here 
Er = h^k1/{2m) is the recoil energy, where kL = 27r/A^ 
for the transverse components and kL = 27r/A|| for the lon- 
gitudinal one. As in the previous section, we use the solution 
of the Mathieu equation to determine the hoppings t^ and t^, 
and we calculate the intraspecies couplings U^^ and U^i in 
the gaussian approximation Ill6ll . We find e.g that the hopping 
ratio t^/t^ spans the interval 0.3 — 0.1 when ]/o,|| goes from 
6 Er to 20 Er', in the same parameter range, the ratios U^^/t^ 
and Uii/ti are above 10, suppressing double occupancy of 
the same species: in this case an appropriate description of the 
system is provided by a model of hardcore bosons, on which 
we will focus our attention later in the paper. The interspecies 
coupling U, which would naturally be very strong in this pa- 
rameter interval (U/t^ > 10), can be suppressed by exploiting 
the above cited Eeshbach resonances. 



-^ = ^^ (1 



27r 



rrii 



Vr 
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(16) 



III. MASS-IMBALANCED TWO-COMPONENT FERMI 
GASES 



These equations indicate the recipe for controlling the interac- 
tions, which can be done by either changing the lattice depth 
with lasers, or by tuning the scattering length using Eeshbach 
resonances. Notice that if the lattice is too deep, multi-band 
processes can occur, but the validity of the one-band approxi- 
mation is generally satisfied in practice. 



A. Overview: population- and mass-imbalanced ID mixtures 

In this section we study the density profiles and correlations 
of a mass- and population imbalanced two-component Eermi 
gas with attractive interactions in a harmonic trap. This ex- 
tends previous studies of population imbalanced mixtures on 



the one hand and mass-imbalanced systems on the other hand. 
We shall briefly review the known results. 

The phase diagram of a spin-imbalanced two-component 
Fermi gas with equal masses and attractive interactions can 
be obtained from the Bethe ansatz, both in the continuum 
1143, |48i] and in lattice models 1461]. Only recently, it was rig- 
orously shown that the partially polarized phase is the one- 
dimensional analogue of the FFLO state by means of exact 



ratio, namely 



2n^ = n; . 



(18) 



numerical methods 11491 IsT 
from mean-field theory ||61 



confirming predictions 

I and bosonization ||63|] . In ID, 
this means that the pair-pair correlation functions are modu- 
lated with Q — kp^ — kpi and decay as a power law, slower 
than any competing correlation in the two-particle channel 
(kpa = Trncr is the Fermi momentum). The computation of 
correlations using the Bethe ansatz is generally very difficult, 
and therefore, the exponents of correlations were first cal- 
culated by numerically solving the corresponding the Bethe 
ansatz equations 11641] . Very recently, first analytical results 
for the FFLO correlations from Bethe ansatz have been pre- 
sented in [65]. In more general models that incorporate the 
coupling of fermions to a molecular channel [|66l - l68l] . the par- 
tially polarized phase is of the FFLO type on the BCS side 
only 163, ItO] • The fate of the FFLO state in coupled one- 
dimensional systems has been recently addressed [|64il7ll - l74|] . 

For spin-imbalanced systems in one dimension, it is by now 
well-established that there is a two- shell structure in a har- 
monic trap 1I47I - I51I] . The central shell is always occupied by 
the partially polarized phase, while in the wings there is a 
fully paired phase at low polarizations and a fully polarized 
phase at large polarizations. Thus, there is a critical polar- 
ization Pc at which the shell structure changes (exactly at Pc, 
the whole system is in a partially polarized phase). This was 
first predicted by applying the local density approximation to 
the exactly known phase diagram ll47i |48|] and then confirmed 
in exact numerical simulations ll49l - [5l|] . The FFLO-type cor- 
relations in the partially polarized phase are stable against the 
presence of a harmonic trap [|51[] . The effect of temperature on 
the shell structure was studied using both Bethe ansatz meth- 
ods combined with LDA 11751] and Quantum Monte Carlo sim- 
ulations iSH UM • In a recent experiment with a 3D array of 
one-dimensional tubes III2I] , the theoretical predictions for the 
shell structure were quantitatively verified 0471 |48l |75|] . For 
more details, see the review [38]. 

A mass-imbalanced system of attractively interacting 
fermions is no longer integrable. The phase diagram for the 
population and mass-imbalanced case was therefore obtained 
by using field theory and DMRG calculations (H, UtI H SI 
ItTI] . In Ref. ll44l] . the main focus was on the partially polarized 
phase which shrinks as the mass-imbalance increases due to 
the instability of a strongly mass imbalanced systems against 
a collapse for P < and phase separation for P > ll43l] . 
The partially polarized phase is, similar to the case of a spin- 
imbalanced mixture, of the ID FFLO type ll44l] . A complete 
phase diagram for the case of a majority of heavy particles 
was presented in Ref. \3oi\ , where in particular, the existence 
of trimer-fluid phases was demonstrated. These trimer phases 
can exist if the densities of the two components obey a certain 



This is a necessary condition only; in Ref. 1361] it was shown 
that the trimer gap goes to zero above a critical density ni^c- 
Moreover, in the trimer-fluid phase, the trimer correlations de- 
cay algebraically [45] , while the (s-wave) pairing correlations 
decay exponentially ll36l] . Note that there are in fact many 
more stable multimer bound states at other rational ratios of 
n^ and n^ ll3^[37l|45l], which we will not study in this work. 
As a consequence of the presence of these trimer phases 
(and phases of liquids of even larger compounds), the par- 
tially polarized phase is separated from the vacuum by an ex- 
tended line in the /i-h phase diagram ll36l] . in contrast to the 
population imbalanced system, where this is a critical point 

Another main difference is the broken particle-hole sym- 
metry of the the mass- and spin-imbalanced mixture with re- 
spect to the spin-imbalanced case only. This has immedi- 
ate consequences for the shell structure in a harmonic trap. 
To illustrate this, let us resort to the local density approxi- 
mation (LDA), valid if the trapping potential is sufficiently 
smooth, as it is the case with parabolic potentials whose har- 
monic oscillator length is much larger than the lattice spacing 
or the interparticle distance. In the case of a single species 
and within LDA, the local properties of the trapped system 
can be quantitatively related to those of a bulk system whose 
chemical potential equals the local chemical potential in the 
trap, namely fii = jj. — Va{i — L/2)'^. In the case of a two- 
component mixture the chemical potential is generally species 
dependent, and one can parametrize the two chemical poten- 
tials as /i^ = /i + h/2 and 11^ = 11 — h/2. The LDA assigns 
to local regions of the trap the behavior of a bulk system with 
local average chemical potential jia- -^ fJ^ai, but with fixed 
chemical potential difference (or magnetic field) h = fi^ — ji^. 

In general, to induce a finite polarization into the system the 
magnetic field h has to overcome (in modulus) either one of 
the following gaps: 

A+ = ^(7V^ + 1,7V;-1)-^(7V^,7V;) 

A- = E{N^ - 1, TV; + 1) - E{N^, TV;) . (19) 

where E{N^^ N^) is the ground state energy of a system with 
a given N^ and TV; in the bulk case. The gaps A± are related 
to the spin (pairing) gap Ag as Ag = A+ + A_ . In other 
words, to obtain P 7^ 0, one needs either h > A+ (P > 0) or 
h < A- (P < 0). In the case of equal masses the polarization 
gaps are both equal to As/ 2, and hence they are known to be 
finite through the Bethe Ansatz solution of the model in the 
case of attractive interactions U < [46] . In particular, as 
depicted in Fig.[5l the Bethe Ansatz solution predicts that the 
As is a decreasing function of /i below half filling. Therefore, 
in this situation, applying the LDA to a trapped system a larger 
spin gap is found in the trap wings where the density is lower: 
h being constant across the trap, this means that a polarization 
P induced by h will appear first at the trap center only, moving 
gradually to the trap wings as h (or P) increases. 
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FIG. 5: (color online) Polarization gaps for the one-dimensional 
Hubbard model for hardcore bosons/fermions with equal popula- 
tions, n is the total particle density. Adapted from Ref. p4l] 



The situation changes radically in the case of mass imbal- 
ance. For a two-component Fermi gas the gaps A± were cal- 
culated via DMRG in Ref. |44] (these results carry over to 
the case of a two-component gas of hard-core bosons to be 
discussed later). As could be anticipated, the two polariza- 
tion gaps become unequal, A+ 7^ A_, and, most importantly, 
they are no longer particle-hole symmetric functions of the 
chemical potential (but the spin gap A5 is). In particular A+ 
becomes an increasing function of the chemical potential / 
the density, which is schematically depicted in Fig. [5] This 
implies that, in the mass imbalanced case, an increasing (pos- 
itive) magnetic field h will first polarize the trap wings, while 
the trap center remains unpolarized and exhibits full pairing. 



B. Predicted shell structure in a trap 

Based on the above considerations and the phase diagram 
of Ref. Il36ll valid for t^/ti = 0.3, one can expect a variety 
of shell structures in a harmonic trap using the local density 
approximation. For P < 0, there will be three cases: (i) a 
partially polarized core plus fully paired wings for < \P\ < 
I Pel I , (ii) an extended region with only a partially polarized 
phase at \Pci\ < P < |Pc2|, and (iii) a partially polarized 
core and fully polarized wings at |P| > \Pc2\- In Ref. |45], it 
was demonstrated that at P = —1/3 there is only a partially 
polarized phase where the commensurability condition for the 
existence of trimers is fulfilled locally at every point in the 
trap, confirming the predictions of Ref. [,36J. In the case of 
P > 0, we expect only one critical polarization Pc3, separat- 
ing a three- shell structure from a two shell structure. At small 
P, the equal-density phase occupies the core of the system, 
followed by an intermediate shell that is partially polarized, 
and a fully polarized region in the wings. At large P > Pc3, 
the core is partially polarized, surrounded by a fully polarized 
shell. One goal of our work is to confirm these predictions by 
using numerically exact DMRG simulations II391 I4Q|1. 

We first compute the density and spin-density profiles at 
zero temperature using DMRG to obtain the state diagram of 
such a system and show that this is consistent with the quali- 



tative expectations drawn from the grand-canonical phase di- 
agram 113^ |43] . We devote particular attention to the stability 
of phases with commensurate densities, i.e., the fully paired 
phase and the trimer phase against the presence of the trap. 
In addition, we calculate the pairing correlations in the pres- 
ence of the trap and show that they are of the FFLO type in 
the partially polarized phase. Next, we investigate the effect 
of varying U on the critical polarization separating different 
shell structures at P > 0. Then, we study the behavior in spin- 
dependent traps, motivated by the discussion from Sec. IIIBI 
Finally, we use quantum Monte Carlo simulations to analyze 
the effect of a finite temperature on the shell structure, fo- 
cussing on the stability of the fully paired phase at P > 0. 



C. Relation to experimental parameters 

Let us first make the connection between the symbols used 
here and the experimental set-up more transparent. We envi- 
sion a ^Li-^^K mixture in an optical lattice. We associate 






(20) 
(21) 



We will work with t^ = 0.3t^, which implies to go to a wave- 
length larger than Am (compare Fig.|4]). We avoid the regime 
of very strong mass imbalance t^ <C t^ since in that regime, 
the system is unstable against phase separation (P > 0) 
or a collapse (P < 0) |43]. Since the bare mass ratio is 
^Li/^K ^ 0.05, it is therefore advantageous to tune the 
system to an intermediate ratio of effective masses. For most 
of the simulations we will work with U = — 4t^ (implying 
U « -13.3t;). 

Furthermore we will consider the ratio of the trapping fre- 
quencies, parameterized by 7^ = l^/V^ as a free parameter, 
and we will first discuss r] = 1. The discussion from Sec. IIIBI 
suggests that for the ^Li-^^K mixture, rj ^ 0.4 and we will 
explore this in our DMRG analysis as well. 



D. Density profiles in a harmonic trap at zero temperature 

We investigate the density profiles at T = using DMRG. 
Our DMRG simulations are done on a chain with L = 140 
sites, and using between m = 400 and m = 800 states. 
We consider variable trapping potentials, confining the atomic 
cloud well within the simulation box. In order to label the 
emergent shell structures, we introduce a set of acronyms: PP 
(partially polarized), ED (equal density), FP (fully polarized), 
TP (trimer phase). Then, we use strings of these acronyms 
that describe the shell structure, going from the core to the 
wings. As an example, the label PP-ED stands for a system 
with a PP phase in the core and an ED phase in the wings. 



1. Equal trapping potentials, Vf = Vi, low densities 

In this section, we consider the case of r] = 1, i.e., equal 
trapping potentials for both components. Typical density and 
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FIG. 6: (Color online) Density profiles for P > 0: (a) At low P, 
there are two shells: ED in the core and PP and FP in the wings, (b) 
At large P, PP sits in the core and FP sits in the wings. Parameters: 
t^ = 0.3tt, V = -4tt, A^ = 50, y = 10"^tt and (a) P = 1/5, 
(b) P = 2/5. (solid lines: (n^), dashed lines: local population 
difference -2(5'f )). 
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FIG. 8: (Color online) P < 0, density profiles that match the trimer 
condition n^. = 2n^ (solid lines: (n^), dashed lines: local population 
difference — 6(5'f )). (a) P = —1/3: the trimer condition is fulfilled 
in the entire trap (first shown in Ref. JH). (b) P = -5/21 > -1/3, 
(c) P = —19/51 < —1/3: PP in the core, the trimer condition is 
fulfilled in a thin outer shell. Parameters: t^ — 0.3t^, U — — 4t^, 
Ar = 51, K = 210-^tt. 



FIG. 7: (Color online) Density profiles for P < 0: (a) At low P, 
there are two shells: PP in the core and FP in the wings, (b) At large 
P, PP sits in the core and FP sits in the wings. Parameters: t^ = 
0.3tt, U = -4tt, A^ = 50, K = 2 • 10"^tt and (a) P = -1/10, 
(b) P = —1/2. (solid lines: {m), dashed lines: local population 
difference -2(5'f )). 



spin density profiles obtained from DMRG simulations are 
shown in Figs. [6] and [7] for P > and P < 0, respectively. 
These results confirm our qualitative expectations discussed 
in Sec. IIII Al namely that at P > 0, there are three shells for 
P < Pc3 (ED-PP-FP) and two at large polarizations P > Pcs 
(PP-FP, compare Fig.©. At P < 0, there are two critical po- 
larizations separating the different shell structures from each 
other: (i) PP-ED for P < P^i, (ii) PP for P^ < P < Pc2, 
and (iii) PP-FP for P > Pc2 (compare Fig. [7]). 

The particular case of Pd < P < Pc2 is addressed in 
Fig. [H Right at P = —1/3 and at sufficiently low density, 
the trimer condition 2n^ = n^ is fulfilled in the entire trapped 
cloud (see Fig. Oa) where we reproduce the results from 
Ref. 145]). In this parameter regime, the results of Refs. [l36ll 
and [45] suggest that the ground state is a trimer fluid. Upon 
deviating only slightly from P = —1/3 (which ensures the 
trimer condition globally), the trimer condition is no longer 
fulfilled locally, i.e. 2{nii) t^ (^«t)- "^^i^ demonstrates that 
in order to stabilize the trimer fluid in a harmonic trap one 
needs to fine-tune the global polarization to P = — 1/3. 



2. Large density regime 

While our main interest is in the low density regime, we 
here also include examples of typical density profiles at large 
densities, i.e., where one of the two components has density 
(^icr) = 1- These are presented in Fig. [9] The first two ex- 
amples are for P < and have (n^i) = 1 in the core of the 
system, followed by a thin partially polarized shell. At small 
P < 0, the outer wings are fully paired [see Fig. Ha)] while at 
large P < 0, they are partially polarized [see Fig. Ob)]. In the 
regime of P > 0, first the light atoms form a band insulator 
(^it) = 1. In the particular example shown in Fig. 0^), the 
core is surrounded by a thin PP phase and a broad fully po- 
larized wing. Qualitatively, one- species band insulators with 
(^icr) = 1 form faster at P < than at P > 0, i.e., by increas- 
ing the total particle number or by making the trap tighter. Not 
surprisingly it is energetically favorable to displace the light 
atoms {(J =t) into the outer regions. 



3. Equal trapping potentials: State diagram 

Our results for the shell structure of a mass- and spin- 
imbalanced system at t^ = 0.3t^ and U = — 4t| are sum- 
marized in the state diagram Fig. (TO] We focus on those re- 
gions of polarization P and effective density p = Ny/V in 
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FIG. 10: (Color online) State diagram for V^ 
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FIG. 9: (Color online) Equal trapping potentials V^ = V^: Emer- 
gence of single species band insulators at large densities. Typical 
density profiles at (a) P = -0.6, (b) P = -0.07, (c) P = 0.5 

(Vi = lO-^tt, U = -4tt, h = 0.3tt). 



which no component has formed a band insulator yet (i.e., 
we restrict ourselves to (n^cr) < 1). Evidently, in most of 
the state diagram the partially polarized phase sits in the core 
with the exception of the region < P < Pc3, where the 
equal density phase occupies the central region of the trap. 
It is interesting to emphasize that this does not happen for a 
purely spin-imbalanced system: there, any arbitrarily small 
polarization pushes the equal density phase to the outer wings 
BtI 1481 poll . In that sense, the mass-imbalance stabilizes the 
fully paired/equal density phase. In particular, if one is inter- 
ested in this phase, it is therefore not necessary to fine-tune 
the global polarization to P = 0. This is a consequence of the 
broken particle-hole symmetry in a mass-imbalanced system 
as discussed in Sec. lIIIAl 

The trimer phase, realized at global polarization P = 
— 1/3, extends up to p ^ 0.22 (triangles up in Fig. [T0|. To 
determine this point we have followed the density in the cen- 
ter of the trap and compared it to the critical density of the 
bulk system beyond which the trimer gap closes [1361] , which 
is at n = 0.5 for U = — 4t^. In Ref. 11451] . it was further shown 
that at sufficiently low density and at P = —1/3, the trimer 
correlations in the trap follow a power law. One possibility to 
observe the formation of trimers in an experiment would be 
to open the trap and to let the gas expand in ID in the optical 
lattice (such an experiment was performed in 2D and 3D op- 
tical lattices with balanced two-component Fermi gases |78]). 
Due to the formation of heavy objects, i.e., the trimers, that 
are protected by the associated excitation gap 11361] . the expan- 
sion should be substantially slower at P = —1/3 compared 
to other polarizations. A similar behavior was seen in numer- 
ical simulations for the expansion of two-component Fermi 
gases with a high density of doublons 173, [Sfl] • Our sugges- 



tion therefore is to monitor the expansion velocity as a func- 
tion of both polarization (at low densities p ^ 0.22) and of 
effective density at P = — 1/3. 

Regarding the critical polarizations Pci (i = 1, 2, 3), we ob- 
serve a weak dependence on effective density. Most notably, 
going to small p further stabilizes the ED-PP-FP regime. This 
is consistent with the usual argument that in one dimension, 
pairing at equal densities is more robust at low densities [1381] . 

It is further instructive to discuss the spatial extent of the 
atomic cloud, and of the clouds of each individual compo- 
nent. This follows the analysis of Ref. [1471] for the one- 
dimensional spin-imbalanced Fermi gas, where it was sug- 
gested that the critical polarization can be read off from the 
polarization dependence of the radii. In fact, a similar analy- 
sis was then used in the experimental work on density; profiles 
of one-dimensional spin-imbalanced Fermi gases 



sity p 



We now take a cut at a fixed p = NVV through the state 
diagram and we estimate the extent Rx (x = n, t, |) as the 
region with a finite particle density, where Rn is the cloud 
radius and R^ are the radii for a =t, i- Typical results are 
shown in Fig. [TTJ Several aspects deserve being mentioned. 
First, the cloud is the smallest in the region Pc2 < P < Pci, 
where only a partially polarized phase is present. Secondly, 
for P < Pc2 (where Pc2 is the critical polarization separating 
PP from PP-FP), one observes R^ < R^, while for Pc2 < 
P < 0, one has R^ = R^ = Rn- At Pd and Pc3, however, 
there are no clear features in the dependence of the radii on 
polarization (we have to keep in mind, though, that the particle 
numbers are fairly small). 

While we have mostly focussed onU = — 4t^, it is also 
interesting to study the effect of U on the stability of the ED 
phase. The dependence of Pcs on U is shown in Fig. [121 In- 
creasing \U\ makes the regime in which the ED phase occu- 
pies the core of the system larger, as expected. 
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4. Equal trapping potentials: FFLO correlations in the PP phase 

We now show that in the PP phase the FFLO correlations 
are robust to the presence of the trap. We focus on the region 
P < 0, since at P > 0, the PP phase sits in the core only for 
P > Pc3. In that regime, however, the spin density is typically 
strongly varying across the trap, which disfavors clear signa- 
tures of the FFLO state. In a harmonic trap, the signatures 
of the FFLO state are the cleanest whenever the polarization 
varies slowly with the local chemical potential ||5ll UM • 

In order to demonstrate the presence of FFLO correlations 
we compute the Fourier transform of the pair-pair correlations 
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Here, and throughout the paper, we adopt the same discretiza- 
tion of the k vectors as in a homogeneous system, namely 




FIG. 13: (Color online) Momentum distribution of pairs in the PP 
phase at P < for V^ = V^ ^ 10~%, U = -At^ (t^ = 0.3t^.). 
We display data for P = -4/15, -2/5, -8/15, -2/3, with the 
peak in n^^^^ shifting from left to right, respectively (we do not see 
a clear signature for < |P| < 1/5). Inset: Position Q of the max- 
imum in n^^^^ vs P (squares) and vs vrrieffPeff (circles, see text in 
Sec. lIIID4t . The dashed line is Q = vrrieffPeff , i.e., Eq. ([24]). 



ncr{k = 0) expresses the number of particles with momentum 
k in the interval [— tt/L, tt/L) for a system of size L. 

Our DMRG results are displayed in Fig. [13] Clearly, in the 
curves with |P| > 0, we observe a maximum in n^^^^ at some 
incommensurate momentum Q, indicative of oscillating pair 
correlations. In a homogeneous system in one dimension and 
in the FFLO state, the modulation Q is given by 



Q = kFt- kpi = 7r{N^ - N^)/L . 



(23) 



i.e., Q oc P. In a harmonic trap, first of all, not all majority 
fermions participate in the FFLO state since all particles in the 
FP regions have to be excluded. Second, the length entering 
in Eq. ([23]) is the one that is actually occupied by the quasi- 
condensate. Therefore, Q is not simply proportional to the 
global polarization. If one accounts for that, following the 
procedure described in Ref . [pll] then one obtains 



Q = TrneffPeff 



(24) 



where neff = A^eff/^eff and Peff = (A^t,eff - ^i,eff)/A^eff. 

Leff is the region occupied by the quasi-condensate, and all 
^cr,eff are obtained by integrating (n^cr) over that region. 
Equation ([24]) should properly describe the scaling of the posi- 
tion Q of the maximum in n^^^^ if the spin-density is approx- 
imately constant in the PP shell. We plot Q vs both the total 
polarization P (circles) and Trneff ^eff (squares) in the inset of 
Fig. [13] Q deviates from Eq. (l24l) by up to 10%. This can 
partially be explained by taking into account the uncertainties 
of determining Q and Leff from the finite- size data, and by 
finite- size effects. 

To complete the discussion, one would need to show that 
the FFLO correlations decay algebraically along the trap (see 
the discussion in Ref. |51] for the spin-imbalanced case). For 
a homogeneous system with mass- and spin-imbalance, the 
momentum distribution functions were discussed in Refs. [|43l 
14411 . while in Refs. |[3^|45|], it was shown that away from P = 
— 1/3, the (s-wave) paring correlations decay algebraically, 
modulated with cos{Qx). 
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5. Unequal trapping potentials, V^ / V^ 

We now turn to the case of unequal trapping potentials, r] ^ 
1. T] > 1 implies that the light atoms are squeezed into the core 
of the trapped gas whereas for 7/ < 1, they are being pushed 
out. 

Zero global polarization - In the case of 77 7^ 1, a partially 
polarized phases can be induced if the trapping potentials are 
sufficiently different form each other (results not shown here). 
However, drastically different trapping potentials are neces- 
sary even at small \U\ ^ t^ and the partially polarized phase 
has a strongly varying spin density. This disfavors this set- 
up as a way of realizing the FFLO state, which leaves clearer 
fingerprints whenever the spin density varies slowly with the 
chemical potential 1 7 1 ] . 

P > 0, stability of the ED phase - Next we address the 
question whether unequal trapping potentials drastically alter 
the state diagram or not. As an example, we consider P = 0.2, 
shown in Figs. [T4ta)-(c) for 77 = 1, 10 and r] = 0.4, respec- 
tively. At 77 = 1, the ED phase occupies the core of the system. 
Only a very large r] achieves enough compression of the light 
particles as to wash out completely the ED phase. At 7/ = 0.4 
on the other hand, the central ED shell grows bigger than for 
T] = 1, as more excess t-particles migrate towards the wings. 

To render this observation more quantitative, we plot the 
critical polarization Pcs as a function of r] in the inset of 
Fig- fT4Tb). Pc3 monotonously decreases with increasing 77, 
with a very weak dependence on r] beyond r] ^ 5. It is im- 
portant to stress that for the experimental parameters sketched 
in Sec. Ill B[ the equal density shell is actually stabilized since 

T] <1. 

The trimer phase, realized at global polarization P = — 1/3 
in the entire trap at sufficiently small densities, is remarkably 
stable against varying r]. This is demonstrated in Fig. \T5\ 
where, starting from the parameters of Fig. |8ja), we vary r] 
from T] = 64 to 77 = 0.4. By compressing the trap for the 
light fermions (rj > 1) the cloud size shrinks drastically, thus 
increasing the density in the center of the trap. Eventually, the 
trimer condition Eq. (ITSl ) is no longer fulfilled along the entire 
trap: the heavy fermions prefer to stay in the wings [see the 
example of 77 = 64 in Fig.[T5ta)]. One also drives the core of 
the system out of the regime in which the trimer gap is finite 
(compare Refs. ||3^|45t]). In the opposite regime, 7/ < 1, the 
cloud expands, yet down to 77 ~ 0.1, we still observe a perfect 
match of {rii) = —6{S^), equivalent to 2{ni^) = (n^^). 
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FIG. 14: (Color online) Unequal trapping potentials V^ ^ V^, P = 
0.2 > 0: Typical density profiles at (a) 7/ = 1, (b) r] = 10, (c) 
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E. Density profiles at finite temperatures 

As seen in the previous sections, the formation of an ED 
dimer liquid in the trap center is very robust to the presence of 
a finite, positive polarization P for a large window of charac- 
teristic densities p. Here we probe the robustness of this phe- 
nomenon to the further effect of finite temperatures, by mak- 
ing use of quantum Monte Carlo (QMC) simulations based on 
the Stochastic Series Expansion (SSE) algorithm [1421] . In our 
present study we use a canonical formulation based on double 
directed-loop updates [81]. Fig. [T6l refers to a path through 



the state diagram of Fig. [TOl This corresponds to an the ex- 
perimentally meaningful situation, in which N^ = 24 is kept 
fixed, while N^ is varied, scanning the P axis at variable p. 
The actual trajectory in the P — p plane is indicated in the 
figure. Here t^ = 0.3t^ and U = — 4t^ as in the previous sec- 
tions, while V^ = Vi = 1.5 ' 10~^t^. As a criterion for dimer 
liquid formation in the trap center we require that the local 
magnetization vanishes over the 10 central sites. This criterion 
is met 8itT ^ (actually T = 4-10~^t^, ensuring the elim- 
ination of thermal effects) over a sizable polarization range 
< P ^ 1/3, as already seen in Fig.[TOl As T is increased 
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FIG. 16: (color online) (a) Finite-temperature state diagram of a mass 
imbalanced mixture with t^ = 0.3t^, U = — 4t^, and V^ = V^ = 
1.5 • 10~^t^. The polarization scan is obtained by fixing N^ — 24 
and varying N^. (b) Trajectory in the P — p plane (compare Fig.fTOl) 
described with the above parameters. 



above zero, the polarization range featuring an unpolarized 
dimer liquid (i.e., the ED phase) in the trap core shrinks, 
but it remains sizable up to temperatures T ^ 0.2 t^/ks, 
proving the robustness of ED pairing to realistic conditions. 
Indeed, even though pairing correlations in one-dimensional 
systems become short-ranged as soon as T becomes finite, 
the finite gaps A± that were defined in Eq. ([T9l ) prevent the 
thermally excited majority (t) particles from flowing from the 
trap wings into the center. The parameter region featuring an 
unpolarized core in Fig. [16] shrinks asymmetrically from the 
side of positive polarizations: indeed at higher P the excess 
t particles sitting in the trap wings have a higher potential 
energy, and therefore they need a smaller thermal energy to 
overcome the gaps A± and to flow into the trap center. 



R Heteronuclear, mass-imbalanced Fermi gases: Summary 

So far we have studied the properties of three pairing states 
of a spin- and population imbalanced Fermi gas in a har- 
monic trap, namely the equal density (or dimer fluid) phase, 
the trimer fluid phase and the FFLO state. 

Our main result for the ED phase is that it occupies the core 
of a trapped system over a wide range of positive polariza- 
tions. Unequal trapping potentials which are typical for a het- 
eronuclear system destabilize the ED phase if the confinement 
for the light particles is much stronger than for the heavy ones, 
whereas in the opposite regime, the ED phase is favored. This 
phase is stable against thermal fluctuations as long as they do 
not overcome the polarization gaps A±. 



IV. MASS-IMBALANCED TWO-COMPONENT ROSE 
GASES 

While fermionic mixtures have been shown to support three 
different types of superfluid phases, FFLO, ED and trimer 
phase, which are stable even in the presence of a trapping 
potential, the question naturally arises whether such pairing 
instabilities are indeed possible when considering bosonic bi- 
nary mixtures. The main difference between fermionic and 
bosonic statistics resides in the fact that bosons can (quasi- 
)condense, resulting in a sharp peak at wavevector k = in 
the momentum distribution 



na{k) 



is. 



ik{ri 



-''^^B,{i;j) 



(25) 



where 



Ba{i;j) = {biX,j) 



(26) 



is the one-body density matrix (OBDM), and bia^b]^ are 
bosonic operators. Strictly speaking, the distinctive feature 
of phases dominated by one-body coherence is the power-law 
decay of the OBDM, Eq. ([26l) - which may or may not result in 
a divergent n{k = 0) peak depending on the exponent of the 
power law. Yet in a trapped system such a power-law behav- 
ior is hard to extract, due to the inhomogeneity imposed by 
the trap; and, more importantly, from an experimental point 
of view the most accessible observable is the momentum dis- 
tribution. Hence in the following we will concentrate on the 
na{k = 0) (or condensate) peak, and in particular on its rel- 
ative changes in height and width as one-body coherence is 
suppressed or enhanced. 

The model Hamiltonian in Eq[T7] has been investigated in 
a series of regimes. If the densities, the masses and the inter- 
action parameters are the same for the two species, a paired 
phase appears above a critical value of the interaction ratio 
\U\/Uaa 1 33, 34, 82], that is, one needs sufficiently strong 
intraspecies repulsion in order to open a spin gap which sup- 
presses one-body coherence. The resulting paired phase has 
an exponentially decaying OBDM, and a power-law decaying 
pair-correlation function. 



,t ;.t 



P{z;j) = {b]^blb,^b,^). 



(27) 



A number of numerical and analytical calculations support 
such a picture, which fully recovers the fermionic case in the 
Tonks-Girardeau (or hardcore) limit U^a ^ oo. In the pres- 
ence of a finite mass imbalance, the ED paired phase may 
evolve into a charge-density-wave phase and then into a crys- 
talline phase for specific filling fractions, as explicitly shown 
in the hardcore case in Ref. [1351] . Finally, a collapsed phase 
appears in the strongly attractive regime \U\/Uaa ^ 1, since 
the intraspecies repulsion is not sufficiently strong to prevent 
very large single- site occupancies. 

On the other hand, in the absence of a trap an arbitrary 
imbalance in the populations will lead to a revival of one- 
body coherence, which becomes algebraically decaying (the 
same property applies to fermions). In the case of bosons. 
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finite-momentum pairing is not observed in the presence of 
population imbalance, so that population imbalanced bosonic 
mixtures with attractive interactions do not feature qualita- 
tively different signatures with respect to mixtures of decou- 
pled species. This issue is quite relevant experimentally, given 
that having perfectly balanced populations is essentially im- 
possible: even in the case of homonuclear mixtures, fluctua- 
tions of order \fN^ are generally expected in the population 
of both species. 



A. Pairing in trapped hardcore bosons with mass and 
population imbalance 



Is then bosonic pairing impossible to observe experimen- 
tally in one-dimensional gases? Luckily the answer is "no", 
if one adds the two further ingredients which are the central 
topic of this paper, namely mass imbalance and a trapping 
potential. In the hard-core case, one can extend the LDA ar- 
gument presented in Sec. IIII Al and conclude that, in presence 
of mass imbalance, a dimer fluid is stable in the center of the 
trap up to a finite polarization, depending on t^/t^^ and U. 

This picture is fully confirmed by a direct calculation of 
trapped hardcore bosons, shown in Fig.[T7] There we simulate 
an experimentally meaningful situation of a trapped hardcore 
boson mixture, in which N^^ = 24 is kept fixed, while N^ is 
varied, thereby varying the polarization P continuously. We 
monitor the evolution of the core density nca (averaged over 
the 10 core sites) revealing that a finite window of positive po- 
larizations exists for which the core remains unpolarized, so 
that nc^ = nci- In this situation ED pairing is robust at the 
trap center, and this has strong signatures in the global coher- 
ence properties of the cloud, captured by the momentum dis- 
tribution. Indeed in the polarization window in which the core 
densities are equal, the peak height n|(/c = 0) is strongly sup- 
pressed, and it hardly changes as P decreases, even though N^ 
is increasing. The peak height n^{k = 0) shows two kinks at 
the two boundaries of the polarization window in question. In 
particular the point P = marks a sharp kink-like minimum 
associated with full pairing of all particles across the trap. 
Similar features are revealed in the full width at half maximum 
(FWHM) of the ncr{k) distributions: the polarization window 
with equal core densities shows a strong enhancement of the 
FWHM for the ^-particles, and its boundaries are marked by 
two kinks in the FWHM for the ^-particles. Hence bosonic 
pairing at the trap center is robust to the presence of a finite 
positive polarization, and it bears strong signatures at the level 
of the momentum distribution. 

A detailed calculation of the polarization gap A+ is still 
missing for the softcore boson case. Nonetheless in the fol- 
lowing we will show that for realistic parameters a softcore 
boson mixture in a trap exhibits a similar behavior as well. 
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FIG. 17: (color online) Polarization dependence of core density and 
momentum distribution in a mixture with N^ = 24 and varying N^. 
The system size L = 140 enters into the definition of na{k — 0). 
Other parameters as in Fig. [16] (a) Core density; (b) height of the 
condensate peak; (c) full width at half maximum (FWHM) of the 
condensate peak (in units of the inverse lattice spacing); (d) trajectory 
in the P — p plane. 



B. Trimer formation in hardcore bosonic mixtures: 
momentum distribution signatures 

In this and all of the following sections, we will focus on 
the case of further population imbalance, and in particular, 
the case P < 0, in which, as discussed earlier, ED pairing 
is always absent in the trap center. In this regime trimer for- 
mation can instead appear when P = —1/3. Such a case of 
strongly imbalanced population has not been extensively stud- 
ied yet; in particular, the emergence of trimer formation has 
up to now been discussed only in systems of fermionic atoms 
or dipolar molecules [36, 37, 83, 84]. The key point here is 
to qualitatively and quantitatively understand the competition 
between possible trimer instabilities and the emergence of ei- 
ther phase separation or single species superfluidity as a func- 
tion of the hopping imbalance, which is expected to play a 
prominent role as in the fermionic case. 

We begin our discussion with the case of hardcore bosons, 
moving then to the case of softcore bosons in the next sec- 
tions. Ref.[3^ has shown that, in the grand-canonical ensem- 
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FIG. 18: (color online) Evolution of core density and momentum 
distribution in a mixture with A^^ = 12 and varying Ni; here V^ = 
Vi = 1.8 • 10~^t^, and the other parameters are as in Fig. [Tt] (a) 
Core density; (b) height of the condensate peak; (c) full width at 
half maximum (FWHM) of the condensate peak; (d) trajectory in the 
P — p plane. 



ing to low densities in the system, for which trimer formation 
is a most robust phenomenon 1.36.1 . We observe that, as ex- 
pected from LDA, the polarization —1/3 is achieved in the 
trap center (nc| = 2nct) ^^^Y when the global polarization 
P is exactly at the same value. When the trimer condition on 
the polarization is satisfied, a trimer liquid appears to form, 
as shown by an anomaly in the one-body coherence proper- 
ties, namely a weak suppression in the height of the conden- 
sate peak n^ik = 0) as well as a slight enhancement of its 
width. But the most notable feature is that the formation of 
trimers at polarization P = —1/3 influences the whole evo- 
lution of the one-body coherence properties for nearby values 
of P. Indeed a sharp kink appears in the height of the conden- 
sate peak n^(/c = 0) atP = —1/3, marking a net change of 
slope in the dependence of this quantity on P. In fact a simi- 
lar kink is observed for P = 0, at which the ED paired phase 
spreads throughout the trap. That kink is associated with the 
fact that ^-particles added to the system to give P < will 
not form dimers with ^-particles, giving rise to a strong en- 
hancement of ^-particle coherence. On the other hand, the 
further enhancement of coherence for P < —1/3 shows that 
for —1/3 < P < the ^-particles added to the system have 
the tendency to form trimer bound states with the t-particles, 
and it is only when the trimer formation has been saturated at 
P = — 1/3 that the one-body coherence grows faster with in- 
creasing 7V|, whence the kink at P = —1/3. In principle the 
extra ^-particles added to give P < —1/3 might form quar- 
tets (3 I's and 1 1), which would lead to an another kink in the 
condensate peak with P becomes lower than —1/2 (saturation 
of quartets), and the same for bound states formed from even 
more particles. On the other hand the binding energy of these 
complexes, if finite, will be increasingly low as the size of the 
complex increases, so that the kinks in the ni{k = 0) curve 
will be increasingly weak. 



ble, a mixture of hardcore bosons or fermions with mass im- 
balance t^/t^^ = 0.3 shows a direct transition from a trimer 
liquid phase with fixed polarization P = —1/3 to a vacuum 
phase as the chemical potential is lowered. Hence, within the 
LDA approximation, a trapped mixture cannot realize a trimer 
liquid phase in the trap center unless its polarization is exactly 
at P = —1/3: indeed, if the polarization differs from —1/3, 
the extra particles of either species cannot be accommodated 
in the trap wings, given that a trimer liquid in the trap center 
can only be flanked by an empty region. 

This condition seems to put a very serious limitation to 
the possibility of observing a trimer liquid phase in a real- 
istic trapped system, given that, as mentioned above, a fine 
tuning in the population of the two species is experimentally 
very hard. Nonetheless, strong signatures of the formation of 
a trimer liquid phase at P = —1/3 can be seen in a broader 
range of polarizations, and specifically in the momentum dis- 
tribution. Fig. [18] shows a polarization scan at very low tem- 
perature in a mixture containing A^^ = 12 particles and a vari- 
able number N^ of ^-particles with t^ = 0.3t^. We chose a 
very weak trapping potential (V^ = V; = 1.8*10~^t|), lead- 



C. Trimer binding energy for bosonic mixtures 

In the present section we focus on softcore bosons, and 
we address the experimentally relevant case of ^^Rb mixtures 
confined in species-dependent optical lattices. As discussed 
in Sec.ini this case is an archetypical case aimed at underpin- 
ning the interplay between bosonic soft-coreness and trimer 
pairing. Such settings provide intermediate hopping asymme- 
try and sufficiently large intraspecies repulsion, as required to 
avoid phase separated regimes for large values of U/t^. We 
will start our discussion by considering the effect of a finite in- 
traspecies interaction in the low-density limit, identifying the 
more favorable parameter setup to observe trimer physics in a 
finite, inhomogeneous system, which will then be investigated 
by means of DMRG and QMC simulations. 

Understanding the three-body pairing mechanism in the 
low-density limit sheds light onto the more complicated 
many-body picture: in fact, a finite trimer binding energy with 
magnitude comparable to the hopping rates tcr is a good start- 
ing point to observe the TP phase at the many-body level, 
as has been noticed in the fermionic case ll36l] . At the few- 
body level, the main difference between the fermionic case 



15 



and the bosonic one stems from the finite intraspecies repul- 
sion, which, in a |[/| :^ t(j perturbation picture, allows for 
additional exchange processes in the bosonic case. At a quali- 
tative level, these terms soften the repulsion between a bound 
t — i pair and a ^-particle, thus potentially increasing the 
trimer binding energy by a factor of order c^ t^/(/7|| ~U). 

In order to get quantitative information, we evaluate the 
trimer binding energy via DMRG simulations, by considering 
the quantity 

Et{L) = El{1, 1) + ^^(0, 1) - 2El(1, 2) (28) 

where El{N^j N^) is the ground state energy of a system of 
length L and population N^^ N^ respectively. The trimer bind- 
ing energy, which corresponds to the trimer gap in the zero- 
density limit 1361], is then defined as Et = limi^^oo Et{L), 
and can be extrapolated using a fourth-order algebraic fit of 
simulation data with L = 20, 40, 60, 80, 100. Typical results 
for the parameter range Vo/Er G [6, 12] are shown in Fig. [19] 
(upper panel), where, as a comparison mark, we also present 
results for the hard-core (HC, Uaa = oo) case. The first 
notable feature is that the trimer binding energy of soft-core 
(SC) bosons is always larger than the HC one, thus confirming 
the raw qualitative picture described above; the difference be- 
tween the two notably shrinks at larger lattice depths, where 
the intraspecies repulsion becomes larger (U^^/t^ :^ 56.9 for 
Vq = 12). Moreover, deeper lattices induce larger mass im- 
balance, which by itself is reflected in a larger binding energy 
with respect to the hopping rate t^. 

In order to establish the optimal experimental setting to ob- 
serve trimer physics, one has to take into consideration that a 
finite temperature T in ultracold gases may indeed prevent any 
relevant observation of paired phases with associated binding 
energy smaller that ksT. In the lower panel of Fig. [T9j we 
show the binding energy in nK units in the same parameter 
regime as in the upper panel: in fact, large values of Vb signif- 
icantly decrease the absolute value of Et, making the effect 
of thermal fluctuations more and more relevant. 

Combining the aforementioned arguments, we conclude 
that the best setting where trimer physics may indeed be ob- 
servable is determined in our case by the condition V^/Er ~ 
6; moreover, since in this case U^i c^ U^^ c^ lOt^, we will 
choose in the following U = — 3t^, as larger interspecies at- 
traction may indeed lead to a collapse in an inhomogeneous 
setup. 



D. Density profiles and correlation functions in a harmonic 
trap at zero temperature 

In the following, we employ DMRG calculations in order 
to elucidate the presence and stability of the trimer phase in 
trapped systems, keeping up to 400 states per block and ap- 
plying up to 10 sweeps, with a truncation error in the final 
DMRG step usually smaller than 10~^. We consider equal 
trapping potentials, V^ = V, starting with a commensurate 
polarization P = —1/3 and then study the effects of incom- 
mensurate polarizations. 
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FIG. 19: Trimer binding energy for ^^Rb mixtures as a function of 
the lattice depth Vb (in units of the recoil energy Er. Note that here, 
the recoil energy is identical for the two components as it only de- 
pends on the bare mass). Panel (a): Et in units of t^. In this case, 
deeper lattices lead to larger mass/ratio, and such a larger Et- Panel 
(b): Et in nK units as considered in Sec. |Ill Deeper lattices strongly 
reduce the temperature scale, as evident from the sharp decrease of 
the binding energy in such scale. For all interaction strengths consid- 
ered, the soft-core case has typical energy scales always larger than 
the hard-core one. 



1. Equal trapping potentials, commensurate polarization 



As we have seen in Sec. IIII D[ fermionic gases with ex- 
actly commensurate polarization can stabilize a trimer phase 
over the entire system, as shown from the density profiles 
in Fig. [SI Since the trimer instability is favored in the low- 
density regime, we focused on population regimes where the 
typical total density in the middle of the trap fluctuates around 
(rii) ^ 0.5. In Fig. [20l a typical density profile is shown 
for iV = 21 particles in a L = 140 site system: as in the 
fermionic case, (n^) and —6{S^) coincide within numerical 
errors all over the system, thus providing a first signature of 
trimer formation. Another independent evidence of trimer for- 
mation is then obtained by looking at the correlation function 
decay in the central part of the system; in particular, in the 
trimer phase, both B{x; L/2) and pairing correlations 



D{x;L/2) = {b\^T/2^iL/2h,L/2+xh,L/2+x) 



(29) 
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are expected to decay exponential with x, whilst trimer corre- 
lations 



T{x-L/2) 
{b\ 



^t 



ut 



'T,L/2n,L/2+l"i,L/2n,L/2+x04.,i./2+xy4-,L/2+l+: 



(30) 



are expected to decay algebraically. In Fig. |2T] the decay of 
correlation functions in the middle of the system for the same 
parameter set of Fig. |20] clearly shows that, while both sin- 
gle particle and dimer correlations decay exponentially, the 
trimer correlations do not, further confirming the stability 
of the trimer phase in the inhomogeneous setup; moreover, 
T(x; L/2) does not show any oscillating behavior with x, in 
accordance with the bosonic nature of the three-body bound 
state I 



2. Incommensurate polarization 



As we have seen in Sec. lIIIDl for the case of fermions and 
in Sec. lIVBl for hardcore bosons, even a very small deviation 
from commensurate densities has a drastic effect on trimer 
phases. The consequences of P ^ —1/3 may be even more 
drastic in the bosonic case, which are not subject to the Pauli 
principle and can thus form higher density regions in the mid- 
dle of the trap. 

We have investigated the same parameter regime as in 
Sec. [nil in the presence of a minimal density imbalance, 
5 = P + 1/3 c^ 0.1, which does also represent the typi- 
cal experimental threshold of population control. Similarly 
to the fermionic case, the trimer liquid is fragile with respect 
to incommensurability. At first, one sees that the condition 
2{n^^i) = (n^^i) is not fulfilled in the middle of the trap: as 
shown in Fig. |22l there is a significant departure from com- 
mensurability even for the smallest density imbalance con- 
sidered. Secondly, correlation functions are also strongly af- 
fected: in particular, single-particle correlations change sig- 
nificantly from exponential to algebraic decay, signally the 
emergence of quasi-condensation. In Fig.[23l we compare the 
superfluid correlation of the heavy species at different density 
imbalance for different values of P: in both incommensurate 
cases, the decay in the center of the trap is algebraic. 
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FIG. 20: Polarization and total density distribution in the P = —1/3 
case for attractive ^^Rb bosonic mixtures. Here, A^ = 21. 
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FIG. 21: Correlation functions in the middle of the trap in the 
P = —1/3 case as a function of the distance with respect to the 
trap center: all correlations decay exponentially except for T, point- 
ing to a dominant trimer instability. Here, Aj{Bj) denotes single 
particle correlation function for the t (i) species. 
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FIG. 22: Polarization and total density distribution in the P = 
—7/15 case for attractive ^^Rb bosonic mixtures. The total number 
of particles in the system is A^ = 21. 



E. Momentum distribution signatures for soft-core bosons: 
finite-temperature effects 

Despite the fragility of trimer formation to slight departures 
from the P = —1/3 condition, a strong signature of this phe- 
nomenon is still observed in the evolution of the one-body 
coherence when P is changed around the value —1/3, as dis- 
cussed in Sec. lIVBl Here we discuss the robustness of this ef- 
fect for finite temperatures, and for the softcore case. Figs.[24l 
and [25] show the height of the condensate peak with a vary- 
ing polarization P for a mixture containing N^ = 12 atoms 
and a variable number N^ of |-atoms. Fig. [24] refers to hard- 
core bosons at increasing temperatures, while Fig. [25] shows 
the case of softcore bosons with the same parameters as for 
^^Rb with Fo ^ ^^r. and with [/ = -3t^. In both cases the 
trapping potential is very weak, Va^ = V^ = 1.8*10~^t^. 

We observe that the distinct signatures of formation of 
dimer and trimer liquids - the kinks in ncr(/c = 0) at P = 
and P = —1/3 respectively - are still observed in the softcore 
case, albeit less clearly than in the hardcore case, particularly 
for what concerns the trimer kink. Both for the softcore and 
the hardcore case the kinks appear to be robust at fairly low 
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FIG. 23: Comparison between superfluid correlations of the heavy 
species at different polarizations. 
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FIG. 25: (color online) Polarization dependence of the condensate 
peak at various temperatures for a trapped mixture of softcore bosons 
with A^^ = 12 and increasing N^. Here Uaa = lOtt and U = — 3t^; 
all other parameters are as in Fig. |24] 
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tions increase the trimer binding energy at fixed mass imbal- 
ance, a very small departure from the commensurate condition 
P = —1/3 drastically changes both correlation functions and 
real space distribution of trapped systems. Nevertheless, suf- 
ficiently low temperatures and densities may indeed allow for 
the observation of such exotic liquids in a trapped gas by look- 
ing at the momentum distribution of the heavier component in 
both hard-core and soft-core mixtures, a sharp signature be- 
ing the strong suppression of zero-momentum contribution as 
a function of P. 



FIG. 24: (color online) Polarization dependence of the conden- 
sate peak at various temperatures for a trapped mixture of hardcore 
bosons with A^^ = 12 and increasing N^. Other parameters are as in 
Fig.[l8] 



temperatures, although the trimer kink is quickly rounded off 
as T reaches values in the range 0.05 — 0.1 t^/Zc^. This shows 
clearly that the observation of trimer formation with the cur- 
rent experimental diagnostics of time-of-flight measurements 
is indeed possible as far as the Hamiltonian parameters are 
concerned {e.g. for ^^Rb mixtures), but it requires extreme 
conditions of very low densities and very low temperatures. 
The two latter requirements are actually compatible, as very 
low temperatures can be in principle achieved by further evap- 
orative cooling or algorithmic cooling [1851] of the atoms al- 
ready loaded in the optical lattice. 



F. Summary: Homonuclear two-component Bose gases 

Binary bosonic mixtures represent a valuable setup where 
both dimer and trimer liquid physics in presence of mass im- 
balance leaves strong signatures in the momentum distribution 
accessible to experiments. In case of hard-core intraspecies 
repulsion, equal density pairing in the central part of the trap 
is relatively stable with respect to thermal effects, and may 
be observed in a wide range of polarizations as long as the 
lighter species density is larger, that is P > 0. Trimer physics 
is instead more delicate: why in general soft-core interac- 



V. DISCUSSION AND SUMMARY 

Ultracold atomic gases represent ideal setups to explore un- 
conventional superfluid states of quantum low dimensional 
systems such as FFLO states and superfluids of composite 
particles. In this work, we have investigated binary bosonic 
and fermionic mixtures in the presence of both finite spin- and 
mass-imbalance, taking into consideration typical experimen- 
tal features such as the inhomogeneity induced by trapping 
potentials and a finite temperature. In the first part, we have 
shown how one can modify the effective mass imbalance by 
considering a proper tuning of the underlying optical lattice 
potential, inducing different tunneling rates for the different 
species. In the second part, taking full advantage of the de- 
tailed microscopic Hamiltonian study of SecJIl we have car- 
ried out a combined DMRG and QMC numerical study of 
both Fermi-Fermi and Bose-Bose mixtures for realistic ex- 
perimental parameters. In the fermionic case, a finite mass- 
imbalance helps in stabilizing exactly paired phases in pres- 
ence of spin-imbalance: in sharp contrast with the equal mass 
case, a majority of light particles stabilizes an equal density 
region in the middle of the trap in a broad regime of total den- 
sities and polarizations. This phase is also robust to finite tem- 
perature effects, as corroborated by QMC simulations. The 
partially polarized phase is of the FFLO type, with a clear sig- 
nature in the momentum distribution of pairs. FFLO should 
be easier to see on the P < side, where the heavy fermions 
are the majority species. We can thus conclude that mass- 
imbalanced fermionic mixtures such as ^Li-^^K represent a 
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valuable setup to observe equal-density pairing and FFLO su- 
perfluidity in ID systems under realistic experimental condi- 
tions. Large mass imbalance is known to lead to a richer show- 
case of superfluid states known as composite liquids: as a case 
study, we have investigated the presence and stability of trimer 
liquids in both bosonic and fermionic mass-imbalanced mix- 
tures. While the typical binding energy of trimer bound states 
is smaller than standard pairs, in the case of exactly commen- 
surate densities trimer liquids are indeed robust in inhomoge- 
neous setups even in the presence of slightly different trapping 
potentials. In the bosonic case, an experimental signature of 
such strongly correlated composite liquids is provided by the 
single particle momentum distribution, which displays a sharp 
kink in its polarization dependence even at finite (albeit small) 
temperature. However, once one departs from P = —1/3, the 
trimer liquid is not stable, as signaled by both the density dis- 
tribution and the correlation function decay in the central part 
of the trap. In summary, the experimental realization of trimer 
liquids in binary, mass-imbalanced mixtures represents a chal- 
lenging task, as both low temperature and densities, combined 



with an accurate control on spin-imbalance, will be necessary. 
Nevertheless, such phases may indeed leave appreciable ex- 
perimental signatures even on easily accessible observables 
such as momentum distribution. 
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